"""
天然气物性计算
"""

# '定义状态方程参数
import math
from yangke.common.config import logger
import numpy as np

# '定义特征参数
# 'MWI-摩尔质量；Ei-能量参数；Ki-体积参数；Gi-定位参数；Qi-四极参数;Fi-高位参数；Si-偶极参数；Wi-组合参数
ei, ki, gi, qi, fi, si, wi = (0, 0, 0, 0, 0, 0, 0)

# '定义二元交互作用参数
eij, uij, kij, gij = (0, 0, 0, 0)
# '***********************************************************************************************************

R, TCM, DCM = (0, 0, 0)
U1, G1, Q1, F1, E1, XIJ, EIJ0, GIJ0, BN, BBN = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
D, T, P, DR = (0, 0, 0, 0)

a = [0, 0.1538326, 1.341953, -2.998583, -0.04831228, 0.3757965, -1.589575, -0.05358847, 0.88659463, -0.71023704,
     -1.471722, 1.32185035, -0.78665925, 0.00000000229129, 0.1576724, -0.4363864, -0.04408159, -0.003433888,
     0.03205905, 0.02487355, 0.07332279, -0.001600573, 0.6424706, -0.4162601,
     -0.06689957, 0.2791795, -0.6966051, -0.002860589, -0.008098836, 3.150547, 0.007224479, -0.7057529, 0.5349792,
     -0.07931491, -1.418465, -5.99905E-17,
     0.1058402, 0.03431729, -0.007022847, 0.02495587, 0.04296818, 0.7465453, -0.2919613, 7.294616, -9.936757,
     -0.005399808, -0.2432567, 0.04987016,
     0.003733797, 1.874951, 0.002168144, -0.6587164, 0.000205518, 0.009776195, -0.02048708, 0.01557322, 0.006862415,
     -0.001226752, 0.002850908]

B = [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3,
     3, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 7, 7, 8, 8, 8, 9, 9]

C = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1,
     1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1]

K = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 2, 2, 2, 4, 4, 0, 0, 2, 2, 2, 4, 4, 4, 4, 0, 1, 1, 2, 2, 3, 3, 4, 4,
     4, 0, 0, 2, 2, 2, 4, 4, 0, 2, 2, 4, 4, 0, 2, 0, 2, 1, 2, 2, 2, 2]

U = [0, 0, 0.5, 1, 3.5, -0.5, 4.5, 0.5, 7.5, 9.5, 6, 12, 12.5, -6, 2, 3, 2, 2, 11, -0.5, 0.5, 0, 4, 6, 21, 23, 22,
     -1, -0.5, 7, -1, 6, 4, 1, 9, -13, 21, 8, -0.5, 0, 2, 7, 9, 22, 23, 1, 9, 3, 8, 23, 1.5, 5, -0.5, 4, 7, 3, 0, 1, 0]

G = [0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0,
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0]

Q = [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
     1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1]

F = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0,
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

S = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

w = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

# 各组分摩尔质量
mw = [16.043, 28.0135, 44.01, 30.07, 44.097, 18.0153, 34.082, 2.0159, 28.01, 31.9988, 58.123, 58.123, 72.15,
      72.15, 86.177, 100.204, 114.231, 128.258, 142.285, 4.0026, 39.948]

ei = [0, 151.3183, 99.73778, 241.9606, 244.1667, 298.1183, 514.0156, 296.355, 26.95794, 105.5348, 122.7667,
      324.0689, 337.6389, 365.5999, 370.6823, 402.636293, 427.72263, 450.325022, 470.840891, 489.558373, 2.610111,
      119.6299]

ki = [0, 0.4619255, 0.4479153, 0.4557489, 0.5279209, 0.583749, 0.3825868, 0.4618263, 0.3514916, 0.4533894,
      0.4186954, 0.6406937, 0.6341423, 0.6738577, 0.6798307, 0.7175118, 0.7525189, 0.784955, 0.8152731, 0.8437826,
      0.3589888, 0.4216551]

gi = [0, 0, 0.027815, 0.189065, 0.0793, 0.141239, 0.3325, 0.0885, 0.034369, 0.038953, 0.021, 0.256692, 0.281835,
      0.332267, 0.366911, 0.289731, 0.337542, 0.383381, 0.427354, 0.469659, 0, 0]

qi = [0, 0, 0, 0.69, 0, 0, 1.06775, 0.633276, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

fi = [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

si = [0, 0, 0, 0, 0, 0, 1.5822, 0.39, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

wi = [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

# '给EIJ赋初值
eij = [
    [1, 0.97164, 0.960644, 1, 0.994635, 0.708218, 0.931484, 1.17052, 0.990126, 1, 1.01953, 0.989844, 1.00235,
     0.999268, 1.107274, 0.88088, 0.880973, 0.881067, 0.881161, 1, 1],
    [0.97164, 1, 1.02274, 0.97012, 0.945939, 0.746954, 0.902271, 1.08632, 1.00571, 1.021, 0.946914, 0.973384, 0.95934,
     0.94552, 1, 1, 1, 1, 1, 1, 1],
    [0.960644, 1.02274, 1, 0.925053, 0.960237, 0.849408, 0.955052, 1.28179, 1.5, 1, 0.906849, 0.897362, 0.726255,
     0.859764, 0.855134, 0.831229, 0.80831, 0.786323, 0.765171, 1, 1],
    [1, 0.97012, 0.925053, 1, 1.02256, 0.693168, 0.946871, 1.16446, 1, 1, 1, 1.01306, 1, 1.00532, 1, 1, 1, 1, 1, 1,
     1],
    [0.994635, 0.945939, 0.960237, 1.02256, 1, 1, 1, 1.034787, 1, 1, 1, 1.0049, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [0.708218, 0.746954, 0.849408, 0.693168, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [0.931484, 0.902271, 0.955052, 0.946871, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1.008692, 1.010126, 1.011501, 1.012821,
     1.014089, 1, 1],
    [1.17052, 1.08632, 1.28179, 1.16446, 1.034787, 1, 1, 1, 1.1, 1, 1.3, 1.3, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [0.990126, 1.00571, 1.5, 1, 1, 1, 1, 1.1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1, 1.021, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1.01953, 0.946914, 0.906849, 1, 1, 1, 1, 1.3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [0.989844, 0.973384, 0.897362, 1.01306, 1.0049, 1, 1, 1.3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1.00235, 0.95934, 0.726255, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [0.999268, 0.94552, 0.859764, 1.00532, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1.107274, 1, 0.855134, 1, 1, 1, 1.008692, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [0.88088, 1, 0.831229, 1, 1, 1, 1.010126, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [0.880973, 1, 0.80831, 1, 1, 1, 1.011501, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [0.881067, 1, 0.786323, 1, 1, 1, 1.012821, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [0.881161, 1, 0.765171, 1, 1, 1, 1.014089, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]]

# '给UIJ赋初值
uij = [
    [1, 0.886106, 0.963827, 1, 0.990877, 1, 0.736833, 1.15639, 1, 1, 1, 0.992291, 1, 1.00367, 1.302576, 1.191904,
     1.205769, 1.219634, 1.233498, 1, 1],
    [0.886106, 1, 0.835058, 0.816431, 0.915502, 1, 0.993476, 0.408838, 1, 1, 1, 0.993556, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [0.963827, 0.835058, 1, 0.96987, 1, 1, 1.04529, 1, 0.9, 1, 1, 1, 1, 1, 1.066638, 1.077634, 1.088178, 1.098291,
     1.108021, 1, 1],
    [1, 0.816431, 0.96987, 1, 1.065173, 1, 0.971926, 1.61666, 1, 1, 1.25, 1.25, 1.25, 1.25, 1, 1, 1, 1, 1, 1, 1],
    [0.990877, 0.915502, 1, 1.065173, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [0.736833, 0.993476, 1.04529, 0.971926, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1.028973, 1.033754, 1.038338, 1.042735,
     1.046966, 1, 1],
    [1.15639, 0.408838, 1, 1.61666, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1, 1, 0.9, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1, 1, 1, 1.25, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [0.992291, 0.993556, 1, 1.25, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1, 1, 1, 1.25, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1.00367, 1, 1, 1.25, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1.302576, 1, 1.066638, 1, 1, 1, 1.028973, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1.191904, 1, 1.077634, 1, 1, 1, 1.033754, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1.205769, 1, 1.088178, 1, 1, 1, 1.038338, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1.219634, 1, 1.098291, 1, 1, 1, 1.042735, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1.233498, 1, 1.108021, 1, 1, 1, 1.046966, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]]

# '给KIJ赋初值
kij = [
    [1, 1.00363, 0.995933, 1, 1.007619, 1, 1.00008, 1.02326, 1, 1, 1, 0.997596, 1, 1.002529, 0.982962, 0.983565,
     0.982707, 0.981849, 0.980991, 1, 1],
    [1.00363, 1, 0.982361, 1.00796, 1, 1, 0.942596, 1.03227, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [0.995933, 0.982361, 1, 1.00851, 1, 1, 1.00779, 1, 1, 1, 1, 1, 1, 1, 0.910183, 0.895362, 0.881152, 0.86752,
     0.854406, 1, 1],
    [1, 1.00796, 1.00851, 1, 0.986893, 1, 0.999969, 1.02034, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1.007619, 1, 1, 0.986893, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1.00008, 0.942596, 1.00779, 0.999969, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0.96813, 0.96287, 0.957828, 0.952441,
     0.948338, 1, 1],
    [1.02326, 1.03227, 1, 1.02034, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [0.997596, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1.002529, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [0.982962, 1, 0.910183, 1, 1, 1, 0.96813, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [0.983565, 1, 0.895362, 1, 1, 1, 0.96287, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [0.982707, 1, 0.881152, 1, 1, 1, 0.957828, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [0.981849, 1, 0.86752, 1, 1, 1, 0.952441, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [0.980991, 1, 0.854406, 1, 1, 1, 0.948338, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]]

# '给GIJ赋初值
gij = [[1, 1, 0.807653, 1, 1, 1, 1, 1.95731, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 0.982746, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [0.807653, 0.982746, 1, 0.370296, 1, 1.67309, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 0.370296, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1.67309, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1.95731, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]]


def transpose(x_i_j):
    """
    二维列表转置

    :param x_i_j:
    :return:
    """
    res = np.array(x_i_j).transpose()
    return res.tolist()


eij = transpose(eij)
uij = transpose(uij)
kij = transpose(kij)
gij = transpose(gij)


def pzo_fdt(D, T, K1, BI, cns):
    global U
    T = T + 273.15
    R = 0.00831451
    DR = D * K1 ** 3
    BMIX = 0
    for n in range(1, 19):
        BMIX = BMIX + BI[n] / T ** U[n]

    Z = 1 + BMIX * D

    for n in range(13, 19, 1):
        Z = Z - DR * cns[n] / T ** U[n]

    for n in range(13, 59, 1):
        Z = Z + cns[n] / T ** U[n] * (B[n] - C[n] * K[n] * DR ** K[n]) * DR ** B[n] * math.exp(-C[n] * DR ** K[n])

    P = D * R * T * Z
    return P, Z, BMIX


class NaturalGas:
    """
    天然气物性计算

    示例：
    natural_gas = NaturalGas(compositions={"CH4": 99, "N2": 1})

    natural_gas.get_density()
    """

    def __init__(self, compositions=None):
        self.compositions = compositions or {}
        self.xx: np.ndarray = self.get_xx()  # 天然气组分体积分数数组，且各组分之和为1
        self.mole_weight = self.get_mole_weight()  # 摩尔质量
        self.xm: np.ndarray = self.get_xm()  # 天然气组分质量分数数组，且各组分之和为1
        self.ratio_CH_number = self.get_ratio_CH()
        self.ratio_CH_mass = self.ratio_CH_number * 12
        self.p = -1  # 压力 MPa
        self.t = -274  # 温度 ℃
        self.z_gas = -1  # 压缩因子
        self.t_refer = -1  # 基准温度 ℃
        self.density_mole = -1  # 摩尔密度，kmol/m3
        self.density = -1  # 密度，kg/m3

    @staticmethod
    def get_components_name():
        return ["CH4", "N2", "CO2", "C2H6", "C3H8", "H2O", "H2S", "H2", "CO", "O2",
                "C4H10", "C4H10-1", "C5H12", "C5H12-1", "C6H14", "C7H16", "C8H18", "C9H20", "C10H22", "He", "Ar"]

    def get_ratio_carbon(self):
        """
        计算单位质量天然气中的碳元素的质量

        :return:
        """
        ratio_comp = [12 / 16, 0, 12 / 44, 24 / 30, 36 / 44, 0, 0, 0, 12 / 28, 0,
                      48 / 58, 48 / 58, 60 / 72, 60 / 72, 72 / 86, 84 / 100, 96 / 114, 108 / 128, 120 / 142, 0, 0]
        return sum(self.xm * np.array(ratio_comp))

    def get_xx(self):
        """
        体积分数数组

        :return:
        """
        x = self.compositions
        xx = [x.get(name, 0) for name in self.get_components_name()]
        # '摩尔分数标准化
        xx = np.array(xx) / sum(xx)
        return xx

    def get_xm(self):
        """
        质量分数数组

        :return:
        """
        xm = self.xx * np.array(mw)
        xm = xm / self.mole_weight
        return np.array(xm)

    def get_mole_weight(self):
        mole_wight = 0  # 摩尔质量，燃气平均分子量
        for i in range(21):
            mole_wight = mole_wight + self.xx[i] * mw[i]
        self.mole_weight = mole_wight
        return self.mole_weight

    def set_pt(self, p, t):
        """
        设置天然气压力和温度

        :param p: MPa
        :param t: ℃
        :return:
        """
        if self.p != p or self.t != t:  # 如果p或t发生变化，则更新
            self.p = p
            self.t = t
        # 更新与p和t有关的参数
        self.z_gas = self.get_z_gas()
        self.density_mole = self.get_density_mole()
        self.density = self.get_density()

    def get_density_mole(self):
        """
        获取摩尔密度
        :return:
        """
        if (self.p == -1) or (self.t < -273.15):
            logger.debug("请先设置燃气温度和压力")
        return self.p / self.get_z_gas(self.p, self.t) / 0.00831451 / (self.t + 273.15)

    def get_density(self):
        """
        获取燃气密度
        :return:
        """
        return self.get_density_mole() * self.mole_weight

    def get_z_gas(self, pressure=None, temperature=None):
        """
        计算天然气压缩因子

        @:param pressure: 压力：MPa
        @:param temperature: 温度，℃
        :return:
        """
        if pressure is not None and temperature is not None:
            if pressure != self.p or temperature != self.t:
                self.set_pt(pressure, temperature)

        xi = [0]
        xi.extend(self.xx)  # 下标从1开始的天然气组分，XI[0]是没有意义的空值
        # '****************************************************计算第二维里系数*************************************************
        k1 = 0
        u1 = 0
        g1 = 0
        q1 = 0
        f1 = 0
        e1 = 0
        for i in range(1, 22):  # 1~21
            k1 = k1 + xi[i] * ki[i] ** 2.5
            u1 = u1 + xi[i] * ei[i] ** 2.5
            g1 = g1 + xi[i] * gi[i]
            q1 = q1 + xi[i] * qi[i]
            f1 = f1 + xi[i] * xi[i] * fi[i]
            e1 = e1 + xi[i] * ei[i]

        # '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
        TCM = 1.261 * e1
        DCM = k1 ** (-1.2)
        k1 = k1 * k1
        u1 = u1 * u1

        for i in range(1, 21):  # 1~20
            for j in range(i + 1, 22):
                xij = xi[i] * xi[j]
                if xij != 0:
                    k1 = k1 + 2 * xij * (kij[i - 1][j - 1] ** 5 - 1) * (ki[i] * ki[j]) ** 2.5
                    u1 = u1 + 2 * xij * (uij[i - 1][j - 1] ** 5 - 1) * (ei[i] * ei[j]) ** 2.5
                    g1 = g1 + xij * (gij[i - 1][j - 1] - 1) * (gi[i] + gi[j])

        k1 = k1 ** 0.2
        u1 = u1 ** 0.2
        # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
        bi = [0]
        for n in range(1, 19):
            bbn = 0
            for i in range(1, 22):
                for j in range(1, 22):
                    xij = xi[i] * xi[j]

                    eij0 = eij[i - 1][j - 1] * (ei[i] * ei[j]) ** 0.5
                    gij0 = gij[i - 1][j - 1] * (gi[i] + gi[j]) / 2

                    BN = (gij0 + 1 - G[n]) ** G[n] * \
                         (qi[i] * qi[j] + 1 - Q[n]) ** Q[n] * \
                         (fi[i] ** 0.5 * fi[j] ** 0.5 + 1 - F[n]) ** F[n] * \
                         (si[i] * si[j] + 1 - S[n]) ** S[n] * \
                         (wi[i] * wi[j] + 1 - w[n]) ** w[n]

                    bbn = bbn + xij * eij0 ** U[n] * (ki[i] * ki[j]) ** 1.5 * BN

            bi.append(a[n] * bbn)

        # '****************************************************求解系数CN********************************************************
        cns = [0] * 59  # 58
        for n in range(13, 59, 1):
            cns[n] = (g1 + 1 - G[n]) ** G[n] * (q1 ** 2 + 1 - Q[n]) ** Q[n] * (f1 + 1 - F[n]) ** F[n] * a[n] * u1 ** U[
                n]

        # '****************************************************开始循环求解**********************************************************
        tol = 0.5 * 10 ** (-9)
        xl = 0.000001
        x2 = 40
        ff1, z, bmix = pzo_fdt(xl, self.t, k1, bi, cns)
        ff2, z, bmix = pzo_fdt(x2, self.t, k1, bi, cns)

        ff1 = ff1 - self.p
        ff2 = ff2 - self.p

        if ff1 * ff2 >= 0:
            return

        # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
        for i in range(5000):
            d = (xl + x2) / 2
            ff, z, bmix = pzo_fdt(d, self.t, k1, bi, cns)
            ff = ff - self.p

            if abs(ff) < tol:
                break
            if ff * ff1 < 0:
                x2 = d
            else:
                xl = d
        self.z_gas = z
        return z

    def get_ncv_mass(self, t_refer=None):
        """
        计算天然气低位发热量，质量发热量
        :param t_refer: 参比温度，℃
        :return: kJ/kg
        """
        allowed_refer_temperature = [0, 15, 20, 25]
        self.t_refer = int(t_refer) if t_refer is not None else self.t_refer
        if self.t_refer not in allowed_refer_temperature:
            logger.debug(f"参比温度只能取值：{allowed_refer_temperature}")
            return -1
        t_refer = self.t_refer
        # 质量发热量
        m_heats = {
            0: np.array([50.043, 0, 0, 47.53, 46.35, 0, 15.19, 119.83, 10.1, 0, 45.57, 45.74, 45.26, 45.36, 45.11,
                         44.93, 44.79, 44.69, 44.61, 0, 0]),
            15: np.array(
                [50.035, 0, 0, 47.52, 46.34, 0, 15.2, 119.91, 10.1, 0, 45.57, 45.72, 45.25, 45.35, 45.11, 44.93,
                 44.79, 44.69, 44.6, 0, 0]),
            20: np.array([50.032, 0, 0, 47.51, 46.34, 0, 9.12, 119.93, 10.1, 0, 45.56, 45.72, 45.05, 45.35, 45.1, 44.92,
                          44.79, 44.69, 44.6, 0, 0]),
            25: np.array([50.029, 0, 0, 47.51, 46.33, 0, 15.2, 119.95, 10.1, 0, 45.56, 45.72, 45.25, 45.35, 45.1, 44.92,
                          44.78, 44.68, 44.6, 0, 0])
        }

        m_heat = m_heats.get(self.t_refer)
        res = sum(self.xm * m_heat) * 1000
        return res

    def get_gcv_mass(self, t_refer=None):
        """
        计算天然气低位发热量，质量发热量，该方法目前支持四个参比温度[0, 15, 20, 25]
        :param t_refer: 参比温度，℃
        :return: kJ/kg
        """
        allowed_refer_temperature = [0, 15, 20, 25]
        self.t_refer = int(t_refer) if t_refer is not None else self.t_refer
        if int(self.t_refer) not in allowed_refer_temperature:
            logger.debug(f"参比温度只能取值：{allowed_refer_temperature}")
            return -1
        # 质量发热量
        m_heats = {
            0: np.array([55.662, 0, 0, 52.02, 50.44, 2.5, 16.52, 142.19, 10.1, 0, 49.45, 49.62, 49.01, 49.1, 48.77,
                         48.53, 48.34, 48.21, 48.09, 0, 0]),
            15: np.array([55.574, 0, 0, 51.95, 50.37, 2.47, 16.5, 141.95, 10.1, 0, 49.39, 49.55, 48.95, 49.04, 48.72,
                          48.47, 48.29, 48.15, 48.04, 0, 0]),
            20: np.array([55.545, 0, 0, 51.93, 50.35, 2.45, 16.5, 141.87, 10.1, 0, 49.37, 49.53, 48.93, 49.03, 48.7,
                          48.45, 48.27, 48.13, 48.02, 0, 0]),
            25: np.array([55.516, 0, 0, 51.9, 50.33, 2.44, 16.49, 141.79, 10.1, 0, 49.35, 49.51, 48.91, 49.01, 48.68,
                          48.44, 48.25, 48.12, 48, 0, 0])
        }
        m_heat = m_heats.get(int(self.t_refer))
        res = sum(self.xm * m_heat) * 1000
        return res

    def get_xh(self, t, t_refer=None):
        """
        计算组分的焓值列表，该方法t_refer可以取任何值

        :param t: ℃
        :param t_refer: ℃
        :return:
        """
        self.t_refer = t_refer if t_refer is not None else self.t_refer
        xh = [self.h_component(name, t, self.t_refer) for name in self.get_components_name()]
        return xh

    def h_component(self, name, t, t_refer=None):
        """
        天然气某组分在温度t下的焓值，参比温度为t_refer，该方法参比温度可以任意取值

        :param name: 组分名
        :param t: ℃
        :param t_refer: ℃
        :return:
        """
        self.t_refer = t_refer if t_refer is not None else self.t_refer
        t_refer = self.t_refer + 273.15
        self.t = t or self.t
        t = self.t + 273.15
        r = 8.31451
        fs = {"CH4": [-176685.0998, 2786.18102, -12.0257785, 0.039176193, -0.0000361905,
                      0.0000000202685, -4.97671E-12, -23313.1436, -74600, 16.0425],
              "C2H6": [-186204.4161, 3406.19186, -19.51705092, 0.075658356, -0.0000820417, 0.0000000506114,
                       -1.31928E-11, -27029.3289, -83851.544, 30.069],
              "C3H8": [-243314.4337, 4656.27081, -29.39466091, 0.118895275, -0.000137631, 0.0000000881482,
                       -2.34299E-11, -35403.3527, -104680, 44.0956],
              "C4H10": [-383446.933, 7000.03964, -44.400269, 0.174618345, -0.00020782, 0.000000133979, -3.55168E-11,
                        -50340.1889, -134990, 58.1222],
              "C4H10-1": [-317587.254, 6176.33182, -38.9156212, 0.158465428, -0.000186005, 0.000000119968,
                          -3.20167E-11, -45403.6339, -125790, 58.1222],
              "C5H12": [-423190.339, 6497.1891, -36.8112697, 0.153242473, -0.000154879,
                        0.000000087499, -2.07055E-11, -51554.1659, -153700, 72.1488],
              "C5H12-1": [-276889.4625, 5834.28347, -36.1754148, 0.153333971, -0.00015284,
                          0.0000000819109, -1.79233E-11, -46653.7525, -146760, 72.1488],
              "C6H14": [-581592.67, 10790.97724, -66.3394703, 0.252371516, -0.000290434, 0.00000018022, -4.61722E-11,
                        -72715.4457, -166920, 86.1754],
              "C7H16": [-581592.67, 10790.97724, -66.3394703, 0.252371516, -0.000290434, 0.00000018022, -4.61722E-11,
                        -72715.4457, -166920, 86.1754],
              "C8H18": [-581592.67, 10790.97724, -66.3394703, 0.252371516, -0.000290434, 0.00000018022, -4.61722E-11,
                        -72715.4457, -166920, 86.1754],
              "C9H20": [-581592.67, 10790.97724, -66.3394703, 0.252371516, -0.000290434, 0.00000018022, -4.61722E-11,
                        -72715.4457, -166920, 86.1754],
              "C10H22": [-581592.67, 10790.97724, -66.3394703, 0.252371516, -0.000290434, 0.00000018022, -4.61722E-11,
                         -72715.4457, -166920, 86.1754],
              "He": [0, 0, 2.5, 0, 0, 0, 0, 745.375, 0, 4.0026],
              "H2": [40783.2321, -800.918604, 8.21470201, -0.012697145, 0.0000175361
                  , -0.0000000120286, 3.36809E-12, 2682.484665, 0, 2.0159],
              "H2S": [9543.80881, -68.7517508, 4.05492196, -0.000301456, 0.0000037685
                  , -0.00000000223936, 3.08686E-13, -3278.45728, -20600, 34.0809],
              "CO": [14890.45326, -292.2285939, 5.72452717, -0.008176235, 0.000014569
                  , -0.0000000108775, 3.02794E-12, -13031.31878, -110535.196, 28.0101],
              "SO2": [-53108.4214, 909.031167, -2.356891244, 0.0220445, -0.0000251078
                  , 0.000000014463, -3.36907E-12, -41137.5212, -296810, 64.0638],
              "Ar": [0, 0, 2.5, 0, 0, 0, 0, -745.375, 0, 39.948],
              "H2O": [-39479.6083, 575.573102, 0.931782653, 0.007222713, -0.00000734256
                  , 0.00000000495504, -1.33693E-12, -33039.7431, -241826, 18.0153],
              "CO2": [49436.5054, -626.411601, 5.30172524, 0.002503814, -0.00000021273
                  , -0.000000000768999, 2.84968E-13, -45281.9846, -393510, 44.0095],
              "O2": [-34255.6342, 484.700097, 1.11901096, 0.004293889, -0.00000068363
                  , -0.00000000202337, 1.03904E-12, -3391.4549, 0, 31.9988],
              "N2": [22103.715, -381.846182, 6.08273836, -0.008530914, 0.0000138465
                  , -0.00000000962579, 2.51971E-12, 710.846086, 0, 28.0134]
              }
        f = fs.get(name)
        h = ((-f[0] / t + f[1] * math.log(t) + f[2] * t + f[3] * t ** 2 / 2 + f[4] * t ** 3 / 3 + f[5] * t ** 4 / 4 + f[
            6] * t ** 5 / 5 + f[7]) * r - f[8]) / f[9] / 2.326
        h_refer = ((-f[0] / t_refer + f[1] * math.log(t_refer) + f[2] * t_refer + f[3] * t_refer ** 2 / 2 + f[
            4] * t_refer ** 3 / 3 + f[
                        5] * t_refer ** 4 / 4 + f[6] * t_refer ** 5 / 5 + f[7]) * r - f[8]) / f[9] / 2.326
        return 2.326 * (h - h_refer)

    def get_h(self, t, t_refer=None):
        """
        燃气显热焓，该方法参比温度可以取任何值
        :param t: 燃气温度，℃
        :param t_refer: 参比温度，℃
        :return:
        """
        self.t_refer = t_refer if t_refer is not None else self.t_refer
        h = self.get_xh(t, t_refer) * self.xm  # 质量焓*质量分数
        return sum(h)

    def get_xianre(self, t, t_refer=None):
        return self.get_h(t, t_refer)

    def get_ratio_CH(self):
        xi = [0]
        xi.extend(self.xx)
        c_number = xi[1] * 1 + xi[2] * 0 + xi[3] * 0 + xi[4] * 2 + xi[5] * 3 + xi[6] * 0 + xi[7] * 0 + xi[8] * 0 + xi[
            9] * 1 + xi[10] * 0 + xi[11] * 4 + xi[12] * 4 + xi[13] * 5 + xi[14] * 5 + xi[15] * 6 + xi[16] * 7 + xi[
                       17] * 8 + xi[18] * 9 + xi[19] * 10 + xi[20] * 0 + xi[21] * 0
        h_number = xi[1] * 4 + xi[2] * 0 + xi[3] * 0 + xi[4] * 6 + xi[5] * 8 + xi[6] * 2 + xi[7] * 2 + xi[8] * 2 + xi[
            9] * 0 + xi[10] * 0 + xi[11] * 10 + xi[12] * 10 + xi[13] * 12 + xi[14] * 12 + xi[15] * 14 + xi[16] * 16 + \
                   xi[17] * 18 + xi[18] * 20 + xi[19] * 22 + xi[20] * 0 + xi[21] * 0
        return c_number / h_number

    def get_ratio_CH_mass(self):
        return self.ratio_CH_number * 12

    def get_ncv_mole(self):
        """
        摩尔低位发热量

        :return: kJ/kmol
        """
        # kJ/kg * g/mol
        return self.get_ncv_mass() * self.mole_weight

    def get_gcv_mole(self):
        """

        :return: kJ/kmol
        """
        return self.get_gcv_mass() * self.mole_weight

    def get_ncv_voln(self, t_combustion=15, t_metering=15):
        """
        标况下体积低位发热量 kJ/kg*kg/Nm3 = kJ/Nm3。以体积计量的发热量与计量温度和压力有很大关系，但国际上使用的压力均已
        统一到0.101325MPa，因此，不在考虑压力对该发热量的影响，但仍需考虑温度对该热值的影响。
        此外，该发热量还与燃烧参比温度有很大关系，也需要考虑燃烧参比温度对该热值的影响。国际上常用的燃烧参比温度和计量参比温度如下：
                            燃烧参比温度℃     计量参比温度
        ISO(国际标准)          15               15
        中国                  20               20
        日本/法国/西班牙        0                0
        加拿大                15               15
        美国/英国             15.56(60℉)      15.56(60℉)
        韩国/澳大利亚          15                0
        俄罗斯                25                20/0
        德国/意大利            25                0

        :return:
        """
        heat_vmole = {
            0: {
                0: np.array(
                    [35.818, 0, 0, 63.76, 91.18, 0, 23.1, 10.777, 12.62, 0, 118.18, 118.61, 145.69, 146, 173.45, 200.87,
                     228.28, 255.74, 283.16, 0, 0]),
            },
            15: {
                0: np.array(
                    [35.812, 0, 0, 63.75, 91.16, 0, 23.11, 10.784, 12.62, 0, 118.16, 118.57, 145.67, 145.98, 173.43,
                     200.84, 228.25, 255.71, 283.13, 0, 0]),

                15: np.array(
                    [33.948, 0, 0, 60.43, 86.42, 0, 21.91, 10.223, 11.96, 0, 112.01, 112.4, 138.09, 138.38, 164.4,
                     19.039, 216.37, 242.4, 268.39, 0, 0]),
            },
            20: {
                20: np.array(
                    [33.367, 0, 0, 59.39, 84.94, 0, 21.53, 10.05, 11.76, 0, 110.09, 110.47, 135.72, 136.01, 161.59,
                     187.13, 212.67, 238.25, 263.8, 0, 0]),
            },
            25: {
                0: np.array(
                    [35.808, 0, 0, 63.74, 91.15, 0, 23.11, 10.788, 12.63, 0, 118.15, 118.56, 145.66, 145.96, 173.41,
                     200.82, 228.23, 255.69, 283.11, 0, 0]),
                20: np.array(
                    [33.365, 0, 0, 59.39, 84.93, 0, 21.53, 10.052, 11.76, 0, 110.08, 110.47, 135.72, 136.01, 161.58,
                     187.12, 212.66, 238.24, 263.79, 0, 0]),
            },
        }
        heat_vol = heat_vmole.get(int(t_combustion)).get(int(t_metering))
        res = sum(self.xx * heat_vol) * 1000
        res = res / self.get_z_gas(0.101325, t_metering)
        return res

    def get_gcv_voln(self, t_combustion=15, t_metering=15):
        """
        高位发热量，其他说明参考get_ncv_voln

        :param t_combustion:
        :param t_metering:
        :return:
        """
        heat_vmole = {
            0: {
                0: np.array(
                    [39.84, 0, 0, 69.79, 99.22, 2.01, 25.12, 12.788, 12.62, 0, 128.23, 128.66, 157.76, 158.07, 187.53,
                     216.96, 246.38, 275.85, 305.29, 0, 0]),
            },
            15: {
                0: np.array(
                    [39.777, 0, 0, 69.69, 99.09, 1.98, 25.09, 12.767, 12.62, 0, 128.07, 128.48, 157.57, 157.87, 187.3,
                     216.7, 246.1, 275.53, 304.94, 0, 0]),
                15: np.array(
                    [37.706, 0, 0, 66.07, 93.94, 1.88, 23.78, 12.102, 11.96, 0, 121.4, 121.79, 149.36, 149.66, 177.55,
                     205.42, 233.28, 261.19, 289.06, 0, 0]),
            },
            20: {
                20: np.array(
                    [37.049, 0, 0, 64.91, 92.29, 1.84, 23.37, 11.889, 11.76, 0, 119.28, 119.66, 146.76, 147.09, 174.46,
                     201.89, 229.22, 256.69, 284.03, 0, 0]),
            },
            25: {
                0: np.array(
                    [39.735, 0, 0, 69.63, 99.01, 1.96, 25.07, 12.725, 12.63, 0, 127.96, 128.73, 157.44, 157.75, 187.16,
                     216.53, 245.91, 275.32, 304.71, 0, 0]),
                20: np.array(
                    [37.029, 0, 0, 64.88, 95.25, 1.83, 23.36, 11.882, 11.76, 0, 119.23, 119.62, 146.7, 146.99, 174.39,
                     201.76, 229.13, 256.59, 283.92, 0, 0]),
            },
        }
        heat_vol = heat_vmole.get(int(t_combustion)).get(int(t_metering))
        res = sum(self.xx * heat_vol) * 1000
        res = res / self.get_z_gas(0.101325, t_metering)
        return res
